import os
import ee
import geopandas as gpd
from lonboard import Map, PolygonLayer, BitmapTileLayer
from gee_redlist.ee_auth import initialize_ee
from gee_redlist.ee_rle import load_yaml, make_eoo, area_km2, get_aoo_grid_projection3 Criterion B Details
This notebook provides step-by-step details on {term}Criterion B calculations.
4 Setup
Import Python modules
Initialize Earth Engine
initialize_ee(project=os.environ['GOOGLE_CLOUD_PROJECT'])5 Analysis
# Define the ecosystem that this notebook is analyzing
ecosystem_code = 'MMR-T1.1.1'# Load the country config
country_config_path = os.environ['PIXI_PROJECT_ROOT'] + '/config/country_config.yaml'
country_config = load_yaml(country_config_path)
# Extract the GEE project path from the country config
gee_project_path = country_config['gee_project_path']
# Extract the class info for the ecosystem
class_info = [x for x in country_config['classified_image']['classes'] if x['code'] == ecosystem_code][0]
print(f'{class_info = }')
ecosystem_image = {
'asset_id': f"{gee_project_path}/{ecosystem_code}/{class_info['id']}",
'pixel_value': class_info['id']
}
print(f'{ecosystem_image = }')class_info = {'id': 52, 'name': 'Tanintharyi island rainforest', 'code': 'MMR-T1.1.1'}
ecosystem_image = {'asset_id': 'projects/goog-rle-assessments/assets/MMR-T1.1.1/52', 'pixel_value': 52}
classified_image_asset_id = f"{country_config['gee_project_path']}/{country_config['classified_image']['asset_id']}"
print(f'{classified_image_asset_id = }')
class_img = (
ee.Image(classified_image_asset_id)
.eq(ecosystem_image['pixel_value'])
.selfMask()
)
print(f'class_img: {class_img.getInfo()}')classified_image_asset_id = 'projects/goog-rle-assessments/assets/mm_ecosys_v7b'
class_img: {'type': 'Image', 'bands': [{'id': 'b1', 'data_type': {'type': 'PixelType', 'precision': 'int', 'min': 0, 'max': 1}, 'dimensions': [14299, 23576], 'crs': 'EPSG:4326', 'crs_transform': [0.0008084837557075694, 0, 89.60991403135986, 0, -0.0008084837557075694, 28.548369897789982]}], 'properties': {'system:footprint': {'type': 'LinearRing', 'coordinates': [[94.66768735253133, 28.548775583027798], [93.22272474215634, 28.548775580809586], [91.41652148658464, 28.548775539452144], [89.60939625454633, 28.54877430628462], [89.60949715757825, 9.485667642440385], [91.41652148658464, 9.486595817478428], [92.86148408075269, 9.486595814990817], [94.30644671824773, 9.486595784089035], [95.39016867380931, 9.48715124036664], [96.47389058147114, 9.486595770655939], [97.91885317497325, 9.486595836813008], [99.36381579356966, 9.486595799973498], [101.1708401374696, 9.48566765739743], [101.17094098108427, 28.548774348455886], [99.00257513499045, 28.548775586792935], [97.55761254104706, 28.54877552452014], [95.75140925342585, 28.548775587887942], [94.66768735253133, 28.548775583027798]]}}}
# Commented out until bug is fixed: https://github.com/developmentseed/lonboard/issues/1064
# # Determine the coordinates for viewing the image
# longitude, latitude = class_img.geometry().centroid().getInfo()['coordinates']
# longitude, latitude6 Extent of occurrence (EOO) (subcriterion B1)
Set the scale (in meters) for reducing the image pixels to polygons. Use the image’s nominal scale unless is is less than 50 meters per pixel.
reduction_scale = max(class_img.projection().nominalScale().getInfo(), 50)
reduction_scale90
Convert the classified image to vectors.
ecosystem_polygons = class_img.updateMask(1).reduceToVectors(
scale=reduction_scale,
geometry=class_img.geometry(),
geometryType='polygon',
maxPixels=1e12,
bestEffort=False
)Create a convex hull that encompasses the ecosystem polygons.
# convexHull() is called twice as a workaround for a bug
# (https://issuetracker.google.com/issues/465490917)
hull = ecosystem_polygons.geometry().convexHull(maxError=1).convexHull(maxError=1)# Note:
# - area() without the projection argument calculates the geodesic area.
# - area(proj=) calculates the planar area. the area is calculated in projected units (10 km x 10 km grid cells)!
# This is different than the javascript version, which calculates the area without reprojection
aoo_area_km2 = hull.area(1, proj=get_aoo_grid_projection()).getInfo() * 1e2
print(f'The area of the hull is {aoo_area_km2:.2f} km²')The area of the hull is 50141.45 km²
aoo_area_km2_geodesic = hull.area().getInfo() * 1e-6
print(f'The geodesic area of the hull is {aoo_area_km2_geodesic:.2f} km²')The geodesic area of the hull is 50336.64 km²
(aoo_area_km2_geodesic - aoo_area_km2) / aoo_area_km20.003892813201885356
7 TRY calculating area in projected coords
# hull_projected = hull.transform(get_aoo_grid_projection(), maxError=1)
# hull_projected.getInfo()# hull_projected.projection().getInfo()# hull_projected.area(1, get_aoo_grid_projection()).getInfo() * 1e28 DEBUG using shapely for convex hull
ecosystem_polygons_gdf = gpd.GeoDataFrame.from_features(
ecosystem_polygons.getInfo()['features'],
crs='EPSG:4326'
)
# Save to file
ecosystem_polygons_gdf.to_parquet('test_ecosystem_polygons.parquet')ecosystem_polygons_gdf| geometry | count | label | |
|---|---|---|---|
| 0 | POLYGON ((97.70122 10.6348, 97.70203 10.6348, ... | 2 | 1 |
| 1 | POLYGON ((97.7182 12.10704, 97.71981 12.10704,... | 25 | 1 |
| 2 | POLYGON ((97.72305 12.11755, 97.72386 12.11755... | 3 | 1 |
| 3 | POLYGON ((97.72467 12.11917, 97.72547 12.11917... | 275 | 1 |
| 4 | POLYGON ((97.72952 12.13372, 97.73194 12.13372... | 140 | 1 |
| ... | ... | ... | ... |
| 600 | POLYGON ((97.67131 12.04722, 97.67211 12.04722... | 2 | 1 |
| 601 | POLYGON ((97.67131 12.03671, 97.67211 12.03671... | 1 | 1 |
| 602 | POLYGON ((97.67211 12.04883, 97.67373 12.04883... | 3 | 1 |
| 603 | POLYGON ((97.67373 12.03509, 97.67616 12.03509... | 146 | 1 |
| 604 | POLYGON ((97.68101 10.46259, 97.68262 10.46259... | 26 | 1 |
605 rows × 3 columns
from shapely.ops import unary_union
hull_geom = unary_union(ecosystem_polygons_gdf.geometry).convex_hull# Reproject the geometry to ESRI:54034 (World Cylindrical Equal Area) and compute area in projected units (square meters)
from pyproj import CRS, Transformer
import shapely
# Define projections
wgs84 = CRS("EPSG:4326")
cea = CRS("ESRI:54034")
project_to_cea = Transformer.from_crs(wgs84, cea, always_xy=True).transform
# Project hull geometry
hull_geom_projected = shapely.ops.transform(project_to_cea, hull_geom)# Area is calculated in square meters (ESRI:54034's units)
hull_geom_projected.area * 1e-650057.86991307719
# Try to project first, then calculate hull
ecosystem_polygons_gdf_projected = ecosystem_polygons_gdf.set_crs(wgs84).to_crs(cea)
hull_geom_projected = ecosystem_polygons_gdf_projected.geometry.union_all().convex_hull
hull_geom_projected.area * 1e-650057.86991307719
diff = (hull_geom_projected.area * 1e-6 - aoo_area_km2 ) / aoo_area_km2
diff-0.0016668938119499528
9 END TEST
9.1 Display EOO Layers
9.2 Ecosystem Tiles
tile_url = class_img.getMapId(
vis_params={
'palette': ['blue'],
'opacity': 0.5
}
)['tile_fetcher'].url_format
tile_layer = BitmapTileLayer(
data=tile_url,
tile_size=256,
max_requests=-1,
min_zoom=0,
max_zoom=19,
)9.3 Ecosystem Polygons Layer
ecosystem_polygons_layer = PolygonLayer.from_geopandas(
ecosystem_polygons_gdf,
get_fill_color=[255, 0, 0, 127],
stroked=True,
get_line_width=2,
get_line_color=[0, 0, 0, 150],
)9.4 Ecosystem Hull Layer
type(hull_geom_projected)shapely.geometry.polygon.Polygon
hull_gdf = gpd.GeoDataFrame.from_features(
ee.FeatureCollection(hull).getInfo(),
crs='EPSG:4326'
)
hull_layer = PolygonLayer.from_geopandas(
hull_gdf,
get_fill_color=[0, 0, 255, 63],
stroked=True,
get_line_width=200,
get_line_color=[0, 0, 0, 150],
)Display the map.
m = Map(
layers=[
ecosystem_polygons_layer,
tile_layer,
hull_layer,
],
)
m9.5 Verify that the step-by-step results are consistent
# Direct call to `make_eoo()`
aoo_area_km2_direct_call = area_km2(make_eoo(class_img)).getInfo()
print(f'EOO area: {aoo_area_km2_direct_call:.0f} km²')EOO area: 50337 km²
# Assert the two values are close to each other.
# assert math.isclose(aoo_area_km2, aoo_area_km2_direct_call, abs_tol=1e-4)10 Area of Occupancy (AOO) (subcriterion B2)
The protocol for this adjustment includes the following steps:
- Intersect AOO grid with the ecosystem’s distribution map.
- Calculate extent of the ecosystem type in each grid cell (‘area’) and sum these areas to obtain the total ecosystem area (‘total area’).
- Arrange grid cells in ascending order based on their area (smaller first).
- Calculate accumulated sum of area per cell (‘cumulative area’).
- Calculate ‘cumulative proportion’ by dividing ‘cumulative area’ by ‘total area’ (cumulative proportion takes values between 0 and 1).
- Calculate AOO by counting the number of cells with a ‘cumulative proportion’ greater than 0.01 (i.e. exclude cells that in combination account for up to 1% of the total mapped extent of the ecosystem type).
10.1 Intersect AOO grid with the ecosystem’s distribution map
Load the AOO grid projection
aoo_grid_proj = get_aoo_grid_projection()
aoo_grid_proj.getInfo(){'type': 'Projection',
'wkt': 'PROJCS["World_Cylindrical_Equal_Area", \n GEOGCS["WGS 84", \n DATUM["WGS_1984", \n SPHEROID["WGS 84", 6378137.0, 298.257223563, AUTHORITY["EPSG","7030"]], \n AUTHORITY["EPSG","6326"]], \n PRIMEM["Greenwich", 0.0], \n UNIT["degree", 0.017453292519943295], \n AXIS["Longitude", EAST], \n AXIS["Latitude", NORTH]], \n PROJECTION["Cylindrical_Equal_Area"], \n PARAMETER["central_meridian", 0.0], \n PARAMETER["standard_parallel_1", 0.0], \n PARAMETER["false_easting", 0.0], \n PARAMETER["false_northing", 0.0], \n UNIT["m", 1.0], \n AXIS["Easting", EAST], \n AXIS["Northing", NORTH], \n AUTHORITY["ESRI","54034"]]',
'transform': [10000, 0, 0, 0, 10000, 0]}
Extract the grid scale parameters
aoo_x_scale, _, _, _, aoo_y_scale, _ = aoo_grid_proj.getInfo()['transform']
print(f'{aoo_x_scale = } meters')
print(f'{aoo_y_scale = } meters')aoo_x_scale = 10000 meters
aoo_y_scale = 10000 meters
- Create an Earth Engine feature collection of AOO grid cells that intersect with the ecosystem, and calculate the fractional coverage of the ecosystem within the grid cell.
fractional_coverage_fc = class_img.unmask().reduceRegions(
collection=class_img.geometry().coveringGrid(aoo_grid_proj),
reducer=ee.Reducer.mean(),
).filter(ee.Filter.gt('mean', 0))
# Convert the Earth Engine feature collection to a GeoPandas GeoDataFrame.
fractional_coverage_gdf = ee.data.computeFeatures({
"expression": fractional_coverage_fc,
"fileFormat": "GEOPANDAS_GEODATAFRAME",
})
fractional_coverage_gdf.rename(columns={"mean": "coverage"}, inplace=True)
# Set the CRS (Earth Engine data is in EPSG:4326)
fractional_coverage_gdf = fractional_coverage_gdf.set_crs('EPSG:4326')aoo_grid_cell_count = len(fractional_coverage_gdf)
aoo_grid_cell_count206
10.2 Calculate grid cell area and the total ecosystem area
- Calculate extent of the ecosystem type in each grid cell (‘area’) and sum these areas to obtain the total ecosystem area (‘total area’).
fractional_coverage_gdf['area'] = fractional_coverage_gdf['coverage'] * aoo_x_scale * aoo_y_scale
fractional_coverage_gdf.sort_values(by="area")[0:4]| geometry | coverage | area | |
|---|---|---|---|
| 175 | POLYGON ((98.09603 12.57845, 98.18586 12.57845... | 0.000035 | 3481.367444 |
| 125 | POLYGON ((98.09603 11.93102, 98.18586 11.93102... | 0.000079 | 7875.324046 |
| 111 | POLYGON ((98.09603 11.74632, 98.18586 11.74632... | 0.000079 | 7880.464175 |
| 74 | POLYGON ((98.54519 11.1008, 98.63502 11.1008, ... | 0.000079 | 7898.091697 |
total_area_km2 = fractional_coverage_gdf['area'].sum() / 1e6
print(f'Total ecosystem area: {total_area_km2:.0f} km²')Total ecosystem area: 1937 km²
10.3 Calculate cumulative area in ordered cells
- Arrange grid cells in ascending order based on their area (smaller first).
- Calculate accumulated sum of area per cell (‘cumulative area’).
fractional_coverage_gdf = fractional_coverage_gdf.sort_values(by="area")
fractional_coverage_gdf["cumulative_area"] = fractional_coverage_gdf["area"].cumsum()
fractional_coverage_gdf.sort_values(by="area").head()| geometry | coverage | area | cumulative_area | |
|---|---|---|---|---|
| 175 | POLYGON ((98.09603 12.57845, 98.18586 12.57845... | 0.000035 | 3481.367444 | 3481.367444 |
| 125 | POLYGON ((98.09603 11.93102, 98.18586 11.93102... | 0.000079 | 7875.324046 | 11356.691491 |
| 111 | POLYGON ((98.09603 11.74632, 98.18586 11.74632... | 0.000079 | 7880.464175 | 19237.155665 |
| 74 | POLYGON ((98.54519 11.1008, 98.63502 11.1008, ... | 0.000079 | 7898.091697 | 27135.247363 |
| 160 | POLYGON ((97.91637 12.39331, 98.0062 12.39331,... | 0.000123 | 12301.457861 | 39436.705224 |
10.4 Calculate the cumulative proportion
- Calculate ‘cumulative proportion’ by dividing ‘cumulative area’ by ‘total area’ (cumulative proportion takes values between 0 and 1).
fractional_coverage_gdf["cumulative_proportion"] = fractional_coverage_gdf["cumulative_area"] / 1e6 / total_area_km2
fractional_coverage_gdf| geometry | coverage | area | cumulative_area | cumulative_proportion | |
|---|---|---|---|---|---|
| 175 | POLYGON ((98.09603 12.57845, 98.18586 12.57845... | 0.000035 | 3.481367e+03 | 3.481367e+03 | 0.000002 |
| 125 | POLYGON ((98.09603 11.93102, 98.18586 11.93102... | 0.000079 | 7.875324e+03 | 1.135669e+04 | 0.000006 |
| 111 | POLYGON ((98.09603 11.74632, 98.18586 11.74632... | 0.000079 | 7.880464e+03 | 1.923716e+04 | 0.000010 |
| 74 | POLYGON ((98.54519 11.1008, 98.63502 11.1008, ... | 0.000079 | 7.898092e+03 | 2.713525e+04 | 0.000014 |
| 160 | POLYGON ((97.91637 12.39331, 98.0062 12.39331,... | 0.000123 | 1.230146e+04 | 3.943671e+04 | 0.000020 |
| ... | ... | ... | ... | ... | ... |
| 108 | POLYGON ((98.45536 11.65401, 98.54519 11.65401... | 0.549146 | 5.491457e+07 | 1.689946e+09 | 0.872610 |
| 102 | POLYGON ((98.45536 11.56174, 98.54519 11.56174... | 0.561263 | 5.612625e+07 | 1.746072e+09 | 0.901591 |
| 164 | POLYGON ((98.36552 12.39331, 98.45536 12.39331... | 0.586257 | 5.862569e+07 | 1.804698e+09 | 0.931863 |
| 99 | POLYGON ((98.18586 11.56174, 98.27569 11.56174... | 0.607082 | 6.070816e+07 | 1.865406e+09 | 0.963210 |
| 93 | POLYGON ((98.18586 11.46949, 98.27569 11.46949... | 0.712497 | 7.124972e+07 | 1.936656e+09 | 1.000000 |
206 rows × 5 columns
10.5 Calculate AOO
- Calculate AOO by counting the number of cells with a ‘cumulative proportion’ greater than 0.01 (i.e. exclude cells that in combination account for up to 1% of the total mapped extent of the ecosystem type).
aoo_grid_cells = fractional_coverage_gdf[fractional_coverage_gdf["cumulative_proportion"] > 0.01]
print(f"AOO (number of cells with cumulative proportion > 0.01): {len(aoo_grid_cells)}")AOO (number of cells with cumulative proportion > 0.01): 139
aoo_grid_cells_dropped = fractional_coverage_gdf[fractional_coverage_gdf["cumulative_proportion"] <= 0.01]10.6 Display the layers
aoo_grid_cells_layer = PolygonLayer.from_geopandas(
aoo_grid_cells,
get_fill_color=[0, 255, 0, 63],
)
aoo_grid_cells_layer_dropped = PolygonLayer.from_geopandas(
aoo_grid_cells_dropped,
get_fill_color=[255, 0, 0, 63],
)
m = Map(
layers=[
tile_layer,
aoo_grid_cells_layer,
aoo_grid_cells_layer_dropped
],
)
m10.7 Verify that the step-by-step results are consistent
# # Direct call to `make_aoo()`
# aoo_cells_direct_call = make_aoo(class_img)
# print(f'AOO: {aoo_grid_cell_count_direct_call} grid cells')# assert aoo_grid_cell_count == aoo_grid_cell_count_direct_call11 Criterion B Summary
from IPython.display import Markdown, display
display(
Markdown(f'AOO and EOO were measured as '
f'{aoo_grid_cell_count} 10 x 10 km grid cells '
f'and {aoo_area_km2:.0f} km², respectively. '
f'See [Criterion B Details](mmr-t1-1-1-crit-b-test) for more information. '
f'There is no evidence that suggests this ecosystem should meet the criteria to '
f'be listed as Near Threatened. The ecosystem is assessed as Least Concern under '
f'Criterion B1 and B2. '
f'**Least Concern.**'
)
)AOO and EOO were measured as 206 10 x 10 km grid cells and 50141 km², respectively. See Criterion B Details for more information. There is no evidence that suggests this ecosystem should meet the criteria to be listed as Near Threatened. The ecosystem is assessed as Least Concern under Criterion B1 and B2. Least Concern.